********************************************************************************
* Modelling of representative dataset from Air Products
* Lei Zhang
* 02/06/2014
********************************************************************************

SETS
    i 	"Set of Air Products' plants"		/P1*P5/
    l  	"Set of Competitor's plants"        /PC1*PC4/
    j   "Set of Demand points(markets)"     /D1*D20/
    k   "Set of products"                   /O1/
    t   "Set of time periods"               /T1*T67/;

ALIAS(t,tt);

*-------------------------------------------------------------------------------

SCALAR M        "Big M"               		    /10000/;
SCALAR CM       "Min amount of expansion" 		/90/;
SCALAR RATE     "Discount rate"			  		/0.05/;
SCALAR TAX      "Tax rate"				  		/0.45/;
SCALAR NEXP     "max number of expansion"       /10/;
SCALAR FRACT    "fractor in disjunction"        /0.95/;

$include tabledata

* TABLE DE(j,k,t)       "Demand of product each market at each time period(100ton/90day)"
* TABLE DIS(i,j)        "Distance between Air Products' plants and markets(mile)"
* TABLE DISC(ii,j)      "Distance between competitor's plants and markets(mile)"
* TABLE PRICE(k,t)      "Price of products at different time periods(MM$/100ton)"
* TABLE FALPHA(i,k,t)   "Fixed cost(MM$)"
* TABLE FBETA(i,k,t)    "Variable cost of Air Products(MM$/100ton)"
* TABLE FGAMMA(i,k,t)   "Expansion cost of Air Products(MM$/9000ton)"
* TABLE FGAMMAH(i,k,t)  "First 100ton higher expansion cost part of Air Products(MM$/9000ton)"
* TABLE FT(i,k,t)       "Distribution cost of Air Products(MM$/100ton/mile)"
* TABLE FCALPHA(l,k,t) 	"Fixed cost of competitor(MM$)"
* TABLE FCBETA(l,k,t)   "Variable cost of competitor(MM$/100ton)"
* TABLE FCGAMMA(l,k,t)	"Expansion cost of competitor(MM$/9000ton)"
* TABLE FCT(l,k,t)      "Distribution cost of competitor(MM$/100ton/mile)"

* MM$/9000ton -> MM$/100ton
FGAMMA(i,k,t) = FGAMMA(i,k,t) / 90;
FGAMMAH(i,k,t) = FGAMMAH(i,k,t) / 90;
FCGAMMA(l,k,t) = FCGAMMA(l,k,t) / 90;

*-------------------------------------------------------------------------------

VARIABLES
    npv             "net present value"
    cost 			"cost of all plants"
    ;

BINARY VARIABLES
    x(i,k,t)        "if expansion happend"
    w(i,k,t)        "if first expansion happend at time t"
    xc(l,k,t)       "if expansion happend on competitor's facility"
    u(k,t)          "disjunction"
    ;

POSITIVE VARIABLES
	rev 			"revenue"
	cost1 			"fixed cost"
	cost2 			"production cost"
	cost3 			"investment cost"
	cost4 			"transportation cost"
	
    y(i,j,k,t)      "Amount of product that Air Products sells to market"
    yc(l,j,k,t)     "Amount of product that competitor sells to market"
    c(i,k,t)        "Capacity of Air products facility"
    dc(i,k,t)       "Capacity expansion of Air products"
    z(i,k,t)        "If first expansion happend (if c > 0)"
    cc(l,k,t)       "Capacity of competitor's facility"
    dcc(l,k,t)      "Capacity expansion of competitor's facility"
    ;

    c.fx('P1','O1','T1') = 72;
    c.l('P1',k,t) = 72;
    c.fx('P2','O1','T1') = 225;
    c.l('P2',k,t) = 225;
    c.fx('P3','O1','T1') = 450;
    c.l('P3',k,t) = 450;
    c.fx('P4','O1','T1') = 0;
    c.l('P4',k,t) = 0;
    c.fx('P5','O1','T1') = 0;
    c.l('P5',k,t) = 0;

    cc.fx('PC1','O1','T1') = 270;
    cc.l('PC1',k,t) = 270;
    cc.fx('PC2','O1','T1') = 450;
    cc.l('PC2',k,t) = 450;
    cc.fx('PC3','O1','T1') = 90;
    cc.l('PC3',k,t) = 90;
    cc.fx('PC4','O1','T1') = 630;
    cc.l('PC4',k,t) = 630;

    w.fx('P1',k,'T1') = 1;
    w.fx('P2',k,'T1') = 1;
    w.fx('P3',k,'T1') = 1;

*-------------------------------------------------------------------------------

EQUATIONS

    obj             "Objective function maximize Air Products NPV"
    objin			"Objective function minimize cost of all plants"
    
    r 				"Revnue"
    c1				"fixed cost"
    c2				"production cost"
    c3				"investment cost"
    c4				"transportation cost"

    up1(i,k,t)      "amount of capacity expansion"
    up2(i,k,t)      "Investment decision in Air Products capacity expansion"
    up3(i,k,t)      "When expansion exist, first expansion exist"
    up4(i,k,t)      "relation between z and w"
    up5(i,k,t)      "z(i,k,t) <= 1"
    up6(i,k)		"expansion count"

    lo1(l,k,t)      "amount of capacity expansion"
    lo2(l,k,t)      "Investment decision in Air Products capacity expansion"
    lo3(j,k,t)      "Demand satisfaction for all markets"
    lo4(i,k,t)      "Capacity and supply of Air Products"
    lo5(l,k,t)    	"Capacity and supply of competitor"

    dj11(k,t)       "disjunction"
    dj12(k,t)       "disjunction"
    dj13(k,t)       "disjunction"
    dj21(k,t)       "disjunction"
    dj22(k,t)       "disjunction"
    dj23(k,t)       "disjunction"
    ;

obj..               npv =e= rev-cost1-cost2-cost3-cost4;
r..					rev =e= sum((i,j,k,t), power(1/(1+RATE), ord(t)-1)*PRICE(k,t)*y(i,j,k,t));
c1..				cost1 =e= sum((i,k,t), power(1/(1+RATE), ord(t)-1)*FALPHA(i,k,t)*z(i,k,t));
c2..				cost2 =e= sum((i,k,t), power(1/(1+RATE), ord(t)-1)*FBETA(i,k,t)*sum(j, y(i,j,k,t)));
c3..				cost3 =e= sum((i,k,t), power(1/(1+RATE), ord(t)-1)*(FGAMMAH(i,k,t)*w(i,k,t) + FGAMMA(i,k,t)*dc(i,k,t)));
c4..				cost4 =e= sum((i,j,k,t), power(1/(1+RATE), ord(t)-1)*FT(i,k,t)*DIS(i,j)*y(i,j,k,t));

objin..				cost =e= cost2 + cost4 
                        + sum((l,j,k,t), power(1/(1+RATE), ord(t)-1)*FCT(l,k,t)*DISC(l,j)*yc(l,j,k,t))
                        + sum((l,k,t), power(1/(1+RATE), ord(t)-1)*FCBETA(l,k,t)*sum(j, yc(l,j,k,t)))
                        + sum((l,k,t), power(1/(1+RATE), ord(t)-1)*FCGAMMA(l,k,t)*dcc(l,k,t));

up1(i,k,t)..        dc(i,k,t) =e= CM * x(i,k,t);
up2(i,k,t)$(ord(t)+1 le card(t))..  
					c(i,k,t+1) =e= c(i,k,t) + dc(i,k,t);
up3(i,k,t)..        x(i,k,t) =l= z(i,k,t);
up4(i,k,t)..        z(i,k,t) =e= sum(tt$(ord(t) ge ord(tt)), w(i,k,tt));
up5(i,k,t)..        z(i,k,t) =l= 1;
up6(i,k)..			sum(t, x(i,k,t)) =l= NEXP;

lo1(l,k,t)..        dcc(l,k,t) =e= CM * xc(l,k,t);
lo2(l,k,t)$(ord(t)+1 le card(t))..
                    cc(l,k,t+1) =e= cc(l,k,t) + dcc(l,k,t);
lo3(j,k,t)..        sum(i, y(i,j,k,t)) + sum(l, yc(l,j,k,t)) =e= DE(j,k,t);
lo4(i,k,t)..        sum(j, y(i,j,k,t)) =l= c(i,k,t);
lo5(l,k,t)..      	sum(j, yc(l,j,k,t)) =l= cc(l,k,t);

dj11(k,t)..         FRACT*(sum(i, c(i,k,t)) + sum(l, cc(l,k,t))) 
                        - sum(j, DE(j,k,t)) =g= -M * u(k,t);
dj12(k,t)..         sum(l, xc(l,k,t)) =g= -M * u(k,t);
dj13(k,t)..         sum(l, xc(l,k,t)) =l= M * u(k,t);
dj21(k,t)..         FRACT*(sum(i, c(i,k,t)) + sum(l, cc(l,k,t))) 
                        - sum(j, DE(j,k,t)) + 1e-6 =l= M * (1-u(k,t));
dj22(k,t)..         sum(l, xc(l,k,t))-1 =g= -M * (1-u(k,t));
dj23(k,t)..         sum(l, xc(l,k,t))-1 =l= M * (1-u(k,t));

*-------------------------------------------------------------------------------

option optcr = 0.1;
option reslim = 1e10;

model leader 
	/obj,r,c1,c2,c3,c4,objin,up1,up2,up3,up4,up5,up6,lo1,lo2,lo3,lo4,lo5,
		dj11,dj12,dj13,dj21,dj22,dj23/;
model follower
	/objin,obj,r,c1,c2,c3,c4,lo1,lo2,lo3,lo4,lo5,dj11,dj12,dj13,dj21,dj22,dj23/;

*-------------------------------------------------------------------------------

scalar UB1, LB1, UB2,LB2;
scalar temp1;
scalar temp2;
scalar count /10/;
scalar npv1, npv2;

solve leader using mip maximizing npv;

UB1 = npv.l;
UB2 = cost.l;

c.fx(i,k,t) = c.l(i,k,t);

solve follower using mip minimizing cost;

LB1 = npv.l;
LB2 = cost.l;

*display UB1, LB1, UB2, LB2;

file results /results.dat/;
put results;

for(temp1 = 1 to count,
	for(temp2 = 1 to count,
		c.up(i,k,t) = +INF;
		c.lo(i,k,t) = 0;
		c.fx('P1','O1','T1') = 72;
		c.l('P1',k,t) = 72;
		c.fx('P2','O1','T1') = 225;
		c.l('P2',k,t) = 225;
		c.fx('P3','O1','T1') = 450;
		c.l('P3',k,t) = 450;
		c.fx('P4','O1','T1') = 0;
		c.l('P4',k,t) = 0;
		c.fx('P5','O1','T1') = 0;
		c.l('P5',k,t) = 0;

		npv.up = LB1 + (UB1-LB1)/count*temp1;
		npv.lo = LB1 + (UB1-LB1)/count*(temp1-1);
		cost.up = LB2 + (UB2-LB2)/count*temp2;
		cost.lo = LB2 + (UB2-LB2)/count*(temp2-1);

		solve leader using mip maximizing npv;

		npv1 = npv.l;

		c.fx(i,k,t) = c.l(i,k,t);
		npv.up = +INF;
		npv.lo = -INF;
		cost.up = +INF;
		cost.lo = -INF;

		solve follower using mip minimizing cost;

		npv2 = npv.l;

		put 'npv1: ', npv1;
		put ' npv2: ', npv2;
		put ' m1s: ', leader.modelstat;
		put ' m2s:', follower.modelstat/;
	);
);

